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ABSTRACT 

Cosmological N-Body simulations are used for a variety of applications. Indeed progress 
in the study of large scale structures and galaxy formation would have been very limited 
without this tool. For nearly twenty years the limitations imposed by computing power forced 
simulators to ignore some of the basic requirements for modeling gravitational instability. 
One of the limitations of most cosmological codes has been the use of a force softening length 
that is much smaller than the typical inter-particle separation. This leads to departures from 
collisionless evolution that is desired in these simulations. We propose a particle based method 
with an adaptive resolution where the force softening length is reduced in high density regions 
while ensuring that it remains well above the local inter-particle separation. The method, 
called the Adaptive TreePM, is based on the TreePM code. We present the mathematical 
model and an implementation of this code, and demonstrate that the results converge over a 
range of options for parameters introduced in generalizing the code from the TreePM code. 
We explicitly demonstrate collisionless evolution in collapse of an oblique plane wave. We 
compare the code with the fixed resolution TreePM code and also an implementation that 
mimics adaptive mesh refinement methods and comment on the agreement, and disagreements 
in the results. We find that in most respects the ATreePM code performs at least as well as the 
fixed resolution TreePM in highly over-dense regions, from clustering and number density of 
haloes, to internal dynamics of haloes. We also show that the adaptive code is faster than the 
corresponding high resolution TreePM code. 

Key words: gravitation, methods; N-Body simulations, cosmology: large scale structure of 
the universe 



1 INTRODUCTION 

Large scale structures traced by galaxies are believed to have 
formed by amplification of small perturbations (Peebles 1980; Pea- 
cock 1999; Padmanabhan 2002; Bemardeau et al. 2002). Galaxies 
are highly over-dense systems, matter density p in galaxies is thou- 
sands of times larger than the average density p in the universe. 
Typical density contrast (5 = p/p — 1) in matter at these scales 
in the early universe was much smaller than unity. Thus the prob- 
lem of galaxy formation and the large scale distribution of galaxies 
requires an understanding of the evolution of density perturbations 
from sinall initial values to the large values we encounter today. 

Initial density perturbations were present at all scales that have 
been observed (Spergel et al. 2007; Percival et al. 2007). The equa- 
tions that describe the evolution of density perturbations in an ex- 
panding universe have been known for several decades (Peebles 
1974) and these are easy to solve when the amplitude of pertur- 
bations is small. Once density contrast at relevant scales becomes 
comparable to unity, perturbations becomes non-linear and cou- 
pling with perturbations at other scales cannot be ignored. The 



equation for evolution of density perturbations cannot be solved 
for generic initial conditions in this regime. N-Body simulations 
(e.g., see Efstathiou et al. (1985); Bertschinger (1998); Bagla & 
Padmanabhan (1997); Bagla (2005); Dolag et al. (2008)) are often 
used to study the evolution in this regime. Alternative approaches 
can be used if one requires only a limited amount of informa- 
tion and in such a case either quasi-linear approximation schemes 
(Bernardeau et al. 2002; Zel'Dovich 1970; Gurbatov, Saichev, & 
Shandarin 1989; Matarrese et al. 1992; Brainerd, Scherrer, & Vil- 
lumsen 1993; Bagla & Padmanabhan 1994; Sahni & Coles 1995; 
Hui & Bertschinger 1996) or scaling relations (Davis & Peebles 
1977; Hamilton et al. 1991; Iain, Mo, & White 1995; Kanekar 
2000; Ma 1998; Nityananda & Padmanabhan 1994; Padmanabhan 
et al. 1996; Peacock & Dodds 1994; Padmanabhan 1996; Peacock 
& Dodds 1996; Smith et al. 2003) suffice. However, even the ap- 
proximation schemes and scaling relations must be compared and 
calibrated with simulations before these can be used with confi- 
dence. 

Last three decades have seen a rapid development of tech- 
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niques and computing power for cosmological simulations and tiie 
results of these simulations have provided valuable insight into the 
study of structure formation. The state of the art simulations used 
less than 10^ particles two decades ago (Efstathiou et al. 1988) 
and if the improvement had been due only to computing power 
then the largest simulation possible today should have been around 
10^ particles, whereas the largest simulations done till date used 
more than 10^° particles (Springel et al. 2005). Evidently, devel- 
opment of new methods and optimizations has also played a sig- 
nificant role in the evolution of simulation studies (Efstathiou et al. 
1985; Barnes & Hut 1986; Greengard & Rokhlin 1987; Bouchet & 
Hemquist 1988; Jemigan & Porter 1989; Hemquist 1990; Makino 
1990, 1991; Hernquist, Bouchet, & Suto 1991; Couchman 1991; 
Ebisuzaki et al. 1993; Theuns 1994; Brieu, Summers, & Ostriker 
1995; Suisalu & Saar 1995; Xu 1995; Dubinski 1996; Kravtsov, 
Klypin, & Khokhlov 1997; Macfarland et al. 1998; Bode, Ostriker, 
& Xu 2000; Brieu & Evrard 2000; Dehnen 2000; Knebe, Green, & 
Binney 2001; Springel, Yoshida, & White 2001; Kawai & Makino 
2001 ; Makino 2002; Dehnen 2002; Bagla 2002; Bagla & Ray 2003; 
Makino et al. 2003; Bode & Ostriker 2003; Ray & Bagla 2004; Du- 
binski et al. 2004; Makino 2004; Springel 2005; Merz, Pen, & Trac 
2005; Yoshikawa & Fukushige 2005; Wadsley, Stadel, & Quinn 
2004; Thacker & Couchman 2006; Khandai & Bagla 2008) Along 
the way, code developers have also successfully met the challenge 
posed by the emergence of distributed parallel programming. 

In modeling gravitational clustering, cosmological N-Body 
codes should ensure the following: 

• The universe does not have a boundary. Therefore cosmo- 
logical simulations need to be run with periodic boundary condi- 
tions^. The simulation volume should be large enough for the ef- 
fects of missing modes to be small, (Bagla & Ray 2005; Bagla & 
Prasad 2006; Power & Knebe 2006; Prasad 2007; Bagla, Prasad, & 
Khandai 2008). 

• The mass of each particle in simulations should be much 
smaller than mass scales of interest in the simulation output. This 
is to ensure adequate mass resolution. 

• Each particle in an N-Body simulation represents a very large 
number of particles/objects in the universe. Thus we must ensure 
that pair-wise interaction of N-Body particles is softened at scales 
comparable with the local inter-particle separation. If this is not 
ensured then the resulting two body collisions introduce errors in 
the resulting distribution of particles, (Splinter et al. 1998; Biimey 
& Knebe 2002). 

In spite of the vast improvement in computing power, simulators 
have often had to compromise on one or more of these points. Of- 
ten, errors also creep in due to the approximate methods used for 
solving for force in simulations. A large fraction of current cosmo- 
logical N-Body codes suffer from collisionality or force biasing. 
Force is biased when softening lengths e are much larger than the 
local mean inter-particle separation, fij . Whereas a complementary 
effect, collisions, occur whenever e -C fij. The reader is referred 
to Dehnen (2001) for a detailed discussion on these two effects. 
Codes that adapt their softening lengths to local densities gener- 
ally are mostly of the adaptive mesh refinement type. These use a 
grid for solving the Poisson equation and often have anisotropics 
in force at the scales comparable to the size of a grid cell. Codes 
which use fixed softening lengths are not entirely coUisionless, and. 



An exception are simulations of large spherical volumes that do not suffer 
significant deformation during the course of evolution. 



in highly over-dense regions biasing of force also exists. The TPM, 
(Xu 1995; Bode, Ostriker, & Xu 2000; Bode & Ostriker 2003) is a 
particle based code with a one step adaptive resolution. However, 
in this case the use of the unmodified PM approach for computing 
the long range force introduces significant errors at scales compa- 
rable to the grid. In this work we describe a code which addresses 
all three issues of force anisotropy, collisionality and force bias- 
ing by employing an adaptive softening formalism with the hybrid 
TreePM code. 

The choice of the optimal softening length has been discussed 
at length (Merritt 1996; Athanassoula et al. 2000; Dehnen 2001; 
Price & Monaghan 2007). These studies were carried out in the 
context of isolated haloes in dynamical equilibrium, e.g., Plummer 
and Hernquist profiles, therefore errors could be defined clearly. It 
is also possible to compute and compare various physical quanti- 
ties with analytical expressions derived from the distribution func- 
tion. Dehnen (2001) derived analytical expressions for errors in the 
context of these profiles. This work suggested that the optimal soft- 
ening length must adapt to the local inter-particle separation as a 
function of space and time. Price & Monaghan (2007) have de- 
veloped an energy-momentum conserving formalism with adaptive 
softening and demonstrated that it was superior to fixed softening. 

The evolution of perturbations at small scales depends 
strongly on the mass and force resolution. High force resolutions 
can lead to better modeling of dense haloes, but gives rise to two 
body collisions in regions where the softening length is smaller 
than the local inter particle separation (Splinter et al. 1998; Binney 
& Knebe 2002; Diemand et al. 2004; Binney 2004; El-Zant 2006; 
Romeo et al. 2008). As all particles in very high density regions go 
through such a phase during evolution, any errors arising due to two 
body collisions can potentially effect the structure of high density 
haloes that form. A high force resolution without a corresponding 
mass resolution can also give misleading results as we cannot probe 
shapes of collapsed objects Kuhlman, Melott, & Shandarin (1996). 
In addition, discreteness and stochasticity also limit our ability to 
measure physical quantities in simulations, and these too need to 
be understood properly (Thiebaut et al. 2008; Romeo et al. 2008). 
In all such cases, the errors in modeling is large at small scales. 
It is important to understand how such errors may spread to larger 
scales and affect physical quantities. 

This work is organized as follows. In §2 we describe the for- 
malism for adaptive softening in a cosmological N-body code. In 
§3 we review the TreePM method, which forms the backbone of 
our gravity solver and describe how adaptive softening is imple- 
mented with it. We briefly discuss the performance characteristics 
of ATreePM in §4. We discuss vahdation of the ATreePM code in 
§5 and we conclude in §6. 



2 ADAPTIVE FORCE SOFTENING 
2.1 Formalism 

The aim of any coUisionless N-body code is to self consistently 
evolve the phase- space distribution function (DP) f{r,v,t) under 
its own gravitational force field F(r, t): 

dtf + i^rf).^r + (^.f).F = (1) 
F{r,t) = -gJ Jdyd\y^f{v',^r,t) (2) 

The approach that one takes is to sample the DP, by N phase-space 
points, {ri,Vi}j^j., at initial time t = U. Liouville's theorem 
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then states that evolving the trajectories of these points to any time 
t>ti will be a representation of the DF at that time. Since the sys- 
tem is a collisionless one, one has to suppress artificial two-body 
collisions arising out of interactions between particles which were 
generated for sampling the density field. One therefore assigns a 
finite size to N-Body particles which ensures softening of force at 
small scales, instead of assuming these to be point particles. The 
density field when sampled by point particles, 

p(r) = / dV/(r, v,t) = Y,mj5l{T - rj) (3) 

i=i 

is now smoothed at small scales if we assign a finite size e to every 
particle: 



N 



(4) 



Where W{u, e) is known as the smoothing kernel and we have as- 
sumed that particles are spherical in shape. Here rrij is the mass of 
the j"* particle. We can now integrate the Poisson equation to ob- 
tain the expression for the kernel for computing force and potential. 
Both the quantities are softened at scales below the softening length 
e. We choose to work with the cubic spline kernel (Monaghan & 
Lattanzio 1985) whose expression is given below. Complete ex- 
pressions for the potential and force are given in the Appendix (See 
Eqn.(Al, A2)). 



l-6(f) +6(f) 

2(1 -f)', 
0, 



< f < 0.5 
0.5 < ^ < 1.0 
1.0 < ^ 



(5) 

In the context of individual softening lengths the density at the lo- 
cation of the i*^ particle is given by: 



1992). It has been shown recently (Price & Monaghan 2007), that 
this term plays an important role in N-Body simulations. We study 
the impact of ignoring this term in Cosmological simulations. 

A remedy for momentum non-conservation is to use a sym- 
metrized softening length 



1. 



(7) 



and plug it into the expression for density to re-derive a sym- 
metrized expression for the potential as 4>{rij,tij). This prescrip- 
tion changes the softening length and hence the neighborlist, which 
one has to recompute. Another disadvantage is that with this pre- 
scription for symmetrization, the number of neighbors is not fixed 
for every particle and hence errors in all smoothed estimates are not 
the same for every particle. An alternate method (Hemquist & Katz 
1989) is to symmetrize the kernel itself. 
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The total potential is thus 
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(8) 



(9) 



We can use this to write a Lagrangian which is manifestly symmet- 
ric and this ensures momentum conservation. Energy is conserved 
only if the Ve term is retained in the EOM. 
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mimj<pij 



(10) 



Pin) =J2mjW{ni,ei) = J2mjW{ni,ei) (6) ^ 

j=i j=i 

where ij indicate the particle indices, rij = Ir, — rj| and ej = 
e(ri). The summation in principle can be extended upto infinity 
if the kernel is infinite in extent (e.g. plunmier, gaussian kernels). 
But since such kernels tend to bias the force (Dehnen 2001; Price 
& Monaghan 2007) we will work only with those kernels which 
have compact support, in particular the cubic spline kernel. Such 
kernels ensure that the force is Newtonian beyond the softening 
scale. Un is the number of the nearest nearbors within ei for the 
particle i. In the discussion that follows we assume that this num- 
ber is fixed for every particle and sets the value of the softening 
length ti. We implicitly assume it in our suimnation. Integrating 
the Poisson equation we obtain the Green's function for the poten- 
tial <j)ij = <j)(rij, ei), where the functional form for cf) is given in 
the Appendix (See Eqn.(Al, A2)). 

With the introduction of individual softening lengths for par- 
ticles, the symmetry of the potential is lost and momentum conser- 
vation is violated. Since e(r) is now a local quantity, the EOM is 
incomplete if one takes the expression of force with the fixed soft- 
ening length € replacing it by a local softening length e(r). Hence 
energy conservation also gets violated. This is because the force 
is derived from the potential and with the introduction of a local 
softening length e(r) in the potential, the gradient must also act on 
e(r) giving us an extra term. Traditionally this gmd-e (Ve) term has 
been ignored since in typical applications these were found to be 
subdominant when compared to the usual force (Gingold & Mon- 
aghan 1982; Evrard 1988; Hemquist & Barnes 1990; Monaghan 



The EOM of motion can be derived with this Lagrangian (the 
reader is referred to Price & Monaghan (2007) for details) 

dvi 
dt 



dWiM) ^ dWiijt^) 



(11) 



Where the first term is the standard Newtonian force term (we refer 
to it as the term). The second is the energy conserving Ve term 
which would be zero for fixed softening. Notice that all terms are 
antisymmetric iai,j and hence the total momentum is conserved. 
Here C and Q are: 

d(t>{rij,ei) 



3 



dei 



dWijjei) 
' dei 



(12) 



(13) 



Expressions for ^ and |^ are given in the Appendix (See 
A3,A4 andA5). As e, oc ^ we have ^ 



- 4^ . The term D, 

■ipi 



term ensures that the EOM is accurate to all orders in e (Springel 
& Hemquist 2002). 



3 THE ADAPTIVE TREEPM METHOD 
3.1 The TreePM Method 

The TreePM algorithm (Bagla 2002; Bagla & Ray 2003) is a hybrid 
N-Body method which improves the accuracy and performance of 
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the Barnes-Hut (BH) Tree method (Barnes & Hut 1986) by com- 
bining it with the PM method (Bagla & Padmanabhan 1997; Merz, 
Pen, & Trac 2005; Klypin & Shandarin 1983; Miller 1983; Bouchet 
& Kandrap 1985; Bouchet, Adam, & Pellat 1985; Hockney & East- 
wood 1988). The TreePM method explicitly breaks the potential 
into a short-range and a longe-range component at a scale rs ■ The 
PM method is used to calculate long-range force and the short- 
range force is computed using the BH Tree method. Use of the BH 
Tree for short-range force calculation enhances the force resolution 
as compared to the PM method. 

The gravitational force is divided into a long range and a short 
range part using partitioning of unity in the Poisson equation. 



Ir 



'fik 



^^nGpk_ 

Ir I sr 

+<fik 
k? 
A;2 



cxp (-fc^rj) 

[l - exp (-fe^rj)] 



(14) 
(15) 



Here ipfJ' and ififl' are the short-range and long-range potentials in 
Fourier space, p is the density, G is the gravitational coupling con- 
stant and rs is the scale at which the splitting of the potential is 
done. It has been shown that this particular split between the long 
range and the short range force is optimal amongst a large class 
of suitable functional forms (Bagla & Ray 2003). The long-range 
force is computed in Fourier space with the PM method and the 
short-range force is computed in real space with the Tree method. 
The short range potential and force in real space are: 



sr/ \ 



Gm(f){r, e)erfc 



2r, 



F"'(r.e) = Gmf(r,e)C — 



(16) 
(17) 
(18) 



Here erfc is the complementary error function, (^(r, e) and f (r, e) 
are the usual potential and force kernels, respectively. C{r/rs) 
modifies the softened Newtonian force kernel f (r, e) to the short- 
range force F°'"(r, e). The expression for </>(r, e) and f(r, e) de- 
pend on the kernel W{r, e) used for smoothing and are given for 
cubic spUne in the Appendix. We find that tabulating erfc(u) and 
C{u) and using interpolation to compute these functions is much 
more effective than calculating them every time. The short range 
force is below 1% of the total force at r > 5rs. The short range 
force is therefore computed within a sphere of radius rcut — 5rs 
using the BH tree method. 

The long range force falls below 1% of the total force for r < 
0.5rs. Thus the choice of softening length does not affect the long 
range force in a significant manner as long as the force softening 
is done with a kernel that has a compact support and the softening 
length is below 0.5rs. 

The BH tree structure is built out of cells and particles. Cells 
may contain smaller cells (sub-cells) within them. Sub-cells can 
have even smaller cells within them, or they can contain a particle. 
In three dimensions, each cubic cell is divided into eight cubic sub- 
cells. Cells, as structures, have attributes like total mass, location 
of center of mass and pointers to sub-cells. Particles, on the other 
hand have the usual attributes: position, velocity and mass. 

Force on a particle is computed by adding contribution of 
other particles or of cells. A cell that is sufficiently far away can 



be considered as a single entity and we can add the force due to the 
total mass contained in the cell from its center of mass. If the cell 
is not sufficiently far away then we must consider its constituents, 
sub-cells and particles. Whether a cell can be accepted as a single 
entity for force calculation is decided by the cell acceptance crite- 
rion (CAC). We compute the ratio of the size of the cell Lceii and 
the distance r from the particle in question to its center of mass and 
compare it with a threshold value 

e = < Oc (19) 
r 

The error in force increases with Oc Poor choice of Oc can lead 
to significant errors (Salmon & Warren 1994). Many different ap- 
proaches have been tried for the CAC in order to minimize error 
as well as CPU time usage (Salmon & Warren 1994; Springel, 
Yoshida, & White 2001). The tree code gains over direct summa- 
tion as the number of contributions to the force is much smaller 
than the number of particles. 

The TreePM method is characterized therefore by three pa- 
rameters, T's^Tcut and 9c . For a discussion on the optimum choice of 
these parameters the reader is referred to Bagla & Ray (2003). For 
all our tests we choose conservative values ra = 1.0, Vcut = 5.2rs 
and dc = 0.3 which give errors below 1% in force. All lengths are 
specified in units of the PM grid. 



3.2 The Adaptive IteePM Method 

We choose to implement adaptive softening with a modified 
TreePM code (Khandai & Bagla 2008), which incorporates Barnes 
optimization using groups (Barnes 1990; Makino 1991; Yoshikawa 
& Fukushige 2005), into the TreePM code. In principle one can also 
incorporate a similar formalism for treecodes (Springel, Yoshida, & 
White 2001), P^M codes (Couchman 1991; Couchman, Thomas, 
& Pearce 1995) and other variants like TPM (Xu 1995), TreePM 
(Bagla 2002; Springel 2005) and GOTPM (Dubinski et al. 2004). 
We shall discuss one advantage of using the TreePM code with the 
group optimizations below. 



3.2.1 Estimating the Softening Length 

Our first task is to get an estimate of the softening length for each 
particle. A natural way to extract a local length scale uses the nu- 
merical value of local density. The local number density is related 
to the softening length as: 



e(r) 



n(r)) 



(20) 



n(r) is the number density at the location of the particle and we as- 
sume all particles have the same mass. Here, n„ is a reference num- 
ber and we take it to be the number of neighbors used for estimation 
of the number density. The above equation is implicit and can be 
solved iteratively, see, e.g., Springel & Hemquist (2002); Springel 
(2005); Price & Monaghan (2007) for details. Price & Monaghan 
(2007) have shown that errors are not very sensitive to the exact 
value of rin. We choose rt„ = 32 in our simulations, and also com- 
ment on variation in results with this choice. 

We are using the formalism developed by Price & Monaghan 
(2007) for achieving an adaptive resolution in gravitational interac- 
tions of particles. As the formalism was developed in the context 
of SPH codes, and some of the quantities required can be com- 
puted naturally using the SPH method, the same was used even 
though the gravitational interaction is completely colUsionless. We 
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follow a similar implementation and use methods commonly used 
in SPH simulations, even though there are no hydrodynamical ef- 
fects present in the gravitational interaction being studied here. For 
an overview of SPH methods, please see Lucy (1977); Gingold 
&Monaghan (1977). The SPH methods assign values for functions 
like the density to particles by averaging over nearest neighbors. 
The list of nearest neighbors is an essential requirement for com- 
puting anything using these methods. AU quantities (p, Ve, Q, C, 
etc.) can be computed at runtime with this neighborlist once we 
have converged to a value for e by solving Eqn.(20). We compute 
the neighborlist using linked lists (Hernquist & Katz 1989). 

We put bounds on the maximum and minimum softening 
lengths, emax and emin- A maximum bound is required so that 
force softening is restricted to the short-range force only, we choose 
tmax = ^ in order to ensure that the long range force is of or- 
der 1% of the total force (or smaller) at scales where force soft- 
ening is important. This ensures that any errors arising from non- 
modification of the long range force are smaller than 1%, if one 
puts a lower threshold on the maximum allowed error then the scale 
(-max has to be lowered correspondingly. Alternatively, one can 
work with a larger and then it becomes possible to allow a larger 
imax- A lower bound is also required for the reason that we do not 
want a few isolated highly over-dense regions to dominate the CPU 
time requirements. The value of the lower bound must correspond 
to densities that are much higher than highest over-densities of in- 
terest in the simulation. We set this lower botmd emin = rs/lOOO, 
which corresponds to p ~ 10^° p. The lower bound however is not 
critical to the structure of the code and may be omitted. 

In our implementation a neighbor search is carried out only 
upto emax- Particles which do not have n„ neighbors within emax 
are assigned e, = emax and the spherical top hat (STH) density 
is assigned with the number of neighbors within emax ■ For these 
particles we assign = 1 and (" = 0, which makes their Ve^ — 0. 
A similar assignment is carried out for particles which have < 
emin- The Ve term is calculated before the short-range force so 
that the individual softening lengths needed for short-range force 
calculation are also assigned in the process. 

3-2.2 Memory Requirements 

Even though we use two separate data structures, namely linked 
lists for Ve and tree for V</) in order to compute the total short- 
range force, additional memory requirements compared to TreePM 
are minimal: we require one additional array for storing the soften- 
ing lengths. The Ve term does not require an additional array since 
it is a coinponent of the short-range force and it can be computed at 
run time. This is because the two data structures are never required 
at the same time. We require specification of the largest force soft- 
ening length in a given cell (See the subsection on the cell accep- 
tance criterion below.). This amounts to a single precision array of 
the same size as the number of cells in the tree. 

An advantage of using an analytical splitting of force, in the 
manner TreePM does, is that computation of short-range force does 
not need global data structures. For example one can geometri- 
cally divide regions into smaller regions and construct local trees 
and linked lists in them (just like one would go about doing it in 
a distributed code) and iterate through these regions for comput- 
ing short-range force instead of constructing one global data struc- 
ture for the entire simulation volume for computing the short-range 
force (Dubinski et al. 2004). This reduces memory usage signifi- 
cantly and the dominant part is taken up by the arrays required for 
computing the long range PM force. 



3.2.3 Timestepping Criterion and Integration 

We have implemented a hierarchical time integrator similar to that 
used in GADGET-2 Springel (2005), in which particle trajectories 
are integrated with individual timesteps and synchronized with the 
largest timestep. As we allow the block time step^ to vary with 
time, we work with the so called KDK approach (Kick-Drift-Kick) 
in which velocities are updated in two half steps whereas posi- 
tion is updated in a full step. It can be shown that with a vari- 
able time step, KDK performs better than DKD (Drift-Kick-Drift) 
(see Springel (2005) for details.). We give separate PM (long range, 
global) kicks and Tree (short range, individual) kicks^. The block 

PM 

timestep At is determined by the particle which has the maxi- 

PM 

mum PM acceleration ttmax'- 

=5t('-^\'' (21) 

\ ^max / 

Here St is the dimensionless accuracy parameter. In our implemen- 
tation of the hierarchy of time steps, the smaller time steps differ by 
an integer power (n) of 2 from the largest, block time step. An ar- 
ray is then used to store the value n which determines the timestep 
of the particle. Individual timesteps At, are first calculated: 

AC = 5t(^)''' (22) 

and then the appropriate hierarchy n is chosen depending on this 
value. Here af^ and e; are the modulus of the individual short- 
range acceleration (sum of the V0 and Ve terms) and the softening 
lengths respectively. TreePM has a similar time-stepping criterion 
with ei replace by e. 

The code drifts all the particles with the smallest timestep to 
the next time, where a force computation is done for particles that 
require an updation of velocity (Kick). However the neighborlist 
and individual softening lengths ei are computed for all particles 
at every small timestep. This is because, even though some parti- 
cles do not require a velocity update, their neighbors might require 
one for which they would contribute through their updated soften- 
ing lengths. The Ve term however can be computed only for those 
particles requiring a velocity update. 

Within a given block time step, the smaller time steps are con- 
stant for a given particle. The time step changes across block time 
step and this brings in inaccuracies in evolution of trajectories. It 
is possible, in principle, to ensure that the second order accuracy is 
maintained here. However, we find that the time steps for particles 
change very slowly and this change does not affect trajectories in a 
significant manner. 

The Courant condition is satisfied for the choice of St we use. 
Indeed, we chose St by requiring that two particles in a highly ec- 
centric orbit around each other maintain the trajectory correctly for 
tens of orbits. 

3.2.4 Cell Acceptance Criterion For ATreePM 

The adaptive force resolution formalism requires us to symmetrize 
force between particles that are separated by a distance smaller than 
the larger of the two softening lengths. Without this, the momen- 
tum conservation cannot be ensured. For pairs of particles separated 
by larger distances, there is no need to explicitly symmetrize force 

Same as the largest time step. 
^ Tree kick includes the contribution of the Ve term. 
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Figure 1. Wall clock timing for TreePM (left) as a function of e and ATreePM (right) with Ve (solid) and without Ve (dashed) as a function of n. 



as there is no dependence on the softening length at these scales . 
Thus the cell acceptance criterion needs to be changed within the 
tree part of the code to ensure that for pairs of particles separated by 
the critical distance, forces are computed through pairwise particle- 
particle (PP) interaction. In our implementation, along with parti- 
cles, cells (and hence groups) are assigned softening lengths corre- 
sponding to the particle contained within the cell that has the largest 
softening length. The CAC, that is modified when we go from the 
TreePM to the modified TreePM (Khandai & Bagla 2008) has to be 
changed again to take into account the largest softening lengths for 
particles within cells and groups whose interaction is to be com- 
puted. 



,) 



(23) 



^ceii 1 ^group ^hc softcmng lengths of the cell and group (in the 
sense explained above) respectively, rg^^ is the distance separating 
the centers of mass of the group and the cell. I/^^,, is the size of the 
cell and L^^^^^ is the distance to the furthest particle from the cen- 
ter of mass of the group. This CAC ensures that the interaction of 
particles separated by less than the softening length is computed in 
a direct pairwise manner and hence can be explicitly symmetrized. 
These are also the pairs for which the Ve term needs to be com- 
puted. 



As an aside we would like to note that the Tree method does not con- 
serve momentum explicitly. This is because the tree traversal approximates 
the force due to pairwise interactions, and in the process the pairwise sym- 
metry is lost. In the modified Tree method (Barnes 1990; Makino 1991; 
Yoshikawa & Fukushige 2005; Khandai & Bagla 2008) explicit pairwise 
force is computed for particles within each group. Particles within a group 
have a common interaction list for force due to particles outside the group. 
Exact pairwise PP force is computed explicitly for intra-group particles and 
we have a pairwise symmetry for this component, but for interaction with 
particles outside the group there is no explicit pairwise symmetry and hence 
no explicit momentum conservation. 



4 PERFORMANCE CHARACTERISTICS 
4.1 Timing 

We now look at the wall clock time as a measure of performance be- 
tween different codes. We studied evolution of a power law model 
with n = —I using 128"^ particles up to the stage where the scale 
of non-linearity is 6 grid lengths. More details of the run are given 
in the section on validation of the ATreePM code. Figure I shows 
wall clock time as a function of softening length e for TreePM (left 
panel) and the wall clock time as a function of n„ for the ATreePM 
codes. 

We can qualitatively understand the slope for the TreePM 
curve by looking at the equivalent time-stepping criterion as 
Eqn.(22) for TreePM. 



Atr = St 



1/2 



(24) 



From here the naive expectation is that the time taken should scale 
as . The slope of the curve is in the range —0.55 to —0.65. The 
reason for this small deviation Ues in our use of a hierarchy of time 
steps, where trajectories of all the particles are not updated at every 
time step. However, positions of particles are updated at every time 
step and this operation as well as those related to creating the tree 
structure at every step add an overhead. This overhead becomes 
more and more important at small e where we have many more 
levels of hierarchy realized in a simulation. This leads to steepening 
of the curve from the simple expectation given above. 

For the ATreePM, we expect softening lengths to be larger for 
larger n„. On the other hand, a larger n„ implies a larger neigh- 
borlist and the time taken for setting up the neighborlist increases. 
The second effect is the dominant one and we see that the time 
taken for Adaptive TreePM increases with n„. We see that time 
taken by both variants of ATreePM is similar for n„ — 32 and 48. 
This indicates that the time taken for calculation of the Ve term is 
negligible. There is a difference between the timing for Un = 16, 
as the code with the Ve does not evolve the system correctly: 
this can be seen in all the indicators like the amplitude of cluster- 
ing, mass function, etc., presented in the next section. We see that 
TreePM with e = 0.025 takes 50% more time than ATreePM with 
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Figure 2. Left: Cumulative distribution of softening lengths at the final epoch of an N-Body simulation (See Figure 8 bottom-left). Right: Average square error 
(ASE) for ATreePM with Ve (soUd line) and ATreePM without Ve (dashed Une) as a function of Un 



n„ = 32, whereas the time taken are comparable with ATreePM 
with n„ = 48. 

Since e is assigned by hand in TreePM, the use of a small 
softening length means a considerable number of particles obtain 
a small timestep even though the acceleration for these particles is 
small and does not vary rapidly with time. In a hierarchical integra- 
tor that we use, this means that force computation is done more 
often within a block timestep. In ATreePM on the other hand a 
small softening is only assigned to particles in over-dense regions. 
This saves time in the under-dense and not so over-dense regions 
and the ATreePM code devotes more time to evolving trajectories 
in highly over-dense regions. This is illustrated in Figure 2 where 
we have plotted the cumulative distribution of softening lengths at 
the final epoch in one of the simulations used here (see lower-left 
panel of Figure 8). The softening length was computed here for 
n„ = 32. We find that only 5% of the particles have a soften- 
ing length smaller than the smallest softening length e = 0.025 
used for the fixed resolution simulations. Even at this epoch where 
highly non-linear clustering is seen, nearly half the particles have a 
softening length corresponding to the maximum value of 0.5. This 
shows how we are able to evolve the system in an ATreePM simu- 
lation with a lower computational cost while resolving highly over- 
dense regions. 



4.2 Errors 

In this section we discuss the dependence of errors on rin. In cos- 
mological simulations it is difficult to define errors when softening 
lengths are varied due to the lack of a reference setup. Nevertheless 
we can choose the optimal softening length such that one mini- 
mizes the globally averaged fluctuation in force as the softening 
lengths are varied. Following Price & Monaghan (2007), we define 
average square error: 

r ^ 

ASE{n„) = jrJ2 - + (25) 



N is the total number of particles. C is a normalization constant 
which is taken as with fmax being the largest value of 

force in either runs. An„ is the change in n„ and we choose this to 
be 8. With other values of An„ the overall behaviour of the error 
plot remains qualitatively the same. Figure 2 shows ASE as func- 
tion of n„ for ATreePM with Ve (solid line) and ATreePM without 
Ve (dashed line). These were computed for the same clustered dis- 
tribution of particles. The qualitative behavior of ASE here is same 
as seen in (Price & Monaghan 2007). At small Un ATreePM with 
Ve has larger errors than ATreePM without Ve. This is also re- 
flected in the poor evolution of density fluctuations with n„ = 16 
for the ATreePM with Ve, discussed in the next section. Both vari- 
ants have a minima, the value of error at the minima being lower for 
ATreePM with the Ve term. Another interesting feature is that the 
region where errors are small is fairly broad for the ATreePM with 
the Ve term. For larger n„ the error increases sharply for ATreePM 
without the Ve term. Thus the optimal configuration is the one with 
the Ve term and 20 < n„ < 32. 



5 VALroATION OF THE ADAPTIVE TREEPM CODE 

Cosmological N-Body simulations lack the equivalent of equilib- 
rium distribution functions for haloes, e.g., Plummer halo that may 
be used to validate a new code. Use of such equilibrium distribu- 
tions allows one to quantify errors in a clean manner. However, 
given that the formalism we use is the same as that presented by 
Price & Monaghan (2007), and that their implementation works 
well with such tests gives us confidence in the formalism. In the 
following discussion, we will test the Adaptive TreePM code in a 
variety of ways and look for numerical convergence. We compare 
the performance of the Adaptive TreePM with the fixed resolution 
TreePM. We also study the role of the Ve term and check if drop- 
ping this term leads to any degradation in performance of the code. 
This last point is important as most AMR codes do not have an 
equivalent term in the equation of motion. 
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Figure 3. Top row: Phase plots of the normal (to one of the planes of collapse) component of the velocity vs normal component of particle displacement for 
TreePM with e = ' ^ ' ^ (columns 1-3 respectively). Bottom row: Transverse component of velocity vs Normal component of velocity along one of the 
planes of collapse for the above runs. The scatter plot shows this for a random subset of particles, whereas the line shows the average value in a few bins. 
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Figure 4. Top row: Phase plots of the normal component of particle displacement vs normal component of the velocity for ATreePM with €i < ^ for 
n„ = 16 without Ve, n„ = 16 with Ve and n„ = 32 with Ve (columns 1-3 respectively). Bottom row: Normal component of velocity vs Transverse 
component of velocity along one of the planes of collapse for the above runs. The scatter plot shows this for a random subset of particles, whereas the line 
shows the average value in a few bins. 
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5.1 Two-body collisions: plane wave collapse 

We start by checking whether the adaptive force softening sup- 
presses two body collisions in the Adaptive TreePM code. We re- 
peat a test recommended by Melott et al. (1997) where they study 
the collapse of an oblique plane wave. In this test, the collapse 
should not lead to any transverse motions if the evolution is col- 
Usionless. Melott et al. (1997) had shown that codes where the 
force softening length is much smaller than the inter-particle sepa- 
ration are coUisional and lead to generation of significant transverse 
motions. Some authors (Heitmann et al. 2005) have stated that the 
failure to ensure planar symmetry may not have a one to one corre- 
spondence with coUisionality. However, the test is nevertheless an 
important one and we present the results with the ATreePM code 
here. 

We use exactly the same initial perturbations as used by Melott 
et al. (1997). The simulations are done with 64'^ particles and a 64^ 
grid. We choose = 1 grid length. Other simulation parameters 
are as described in the subsection §3.1 on the TreePM code. The 
output is studied at a = 3, following Melott et al. (1997). 

We first conduct the test with the fixed resolution TreePM 
code. Top row of the Figure 3 shows the phase portrait along the 
direction of collapse for different choices of the force softening 
length. The softening length varies between 0.1 < e < 0.5 in units 
of the mean inter-particle separation. We see that the phase por- 
trait in the multi-stream region is heavily distorted for the smallest 
force softening length but is correct for the largest force softening 
length used here. This reinforces the conclusions of Melott et al. 
(1997) that using a force softening length that is much smaller that 
the mean inter-particle separation leads to two body collisions. This 
point is presented again in the lower row of Figure 3 where we show 
a scatter plot of modulus of the transverse (to the direction of col- 
lapse) velocities and radial velocities. This is shown for a random 
sub-set of all particles. We see that transverse motions are signif- 
icant for the simulation with the smallest force softening length, 
but are under control for the largest force softening length. Thick 
line in the lower panels connects the average modulus of the trans- 
verse velocities in bins of magnitude of radial velocity. The visual 
impression gathered from the scatter plot is reinforced in that with 
decreasing force softening length, we get larger transverse motions. 

Figure 4 presents results of simulations with the same initial 
conditions carried out with the Adaptive TreePM code. We show 
the results for Adaptive TreePM without the Ve term, n„ — 16 
(left colunm); with the Ve term, n„ = 16 (middle column), and, 
with the Ve term, n„ = 32 (middle column). In general we do 
not recommend use of ra„ — 16 due to reasons discussed in the 
preceeding subsection on errors, discussion in the following sub- 
sections, and in Price & Monaghan (2007), but we nevertheless use 
it in order to look for early signs of two body collisions in a simula- 
tion with a relatively small number of particles. The phase portrait 
for all the three Adaptive TreePM runs is a faithful representation of 
the expectations. The transverse motions are suppressed strongly, 
almost to the same level as the TreePM simulation with a force 
softening of e = 0.5. One of the reasons for this is that the highest 
overdensities reached in this experiment do not lead to a consid- 
erable reduction in the force softening length. The other reason is 
that the adaptive nature of the code leads to presence of a sufficient 
number of neighbours within a force softening length and hence the 
anisotropy in the transverse force is reduced substantially. 

Thick lines in the lower panels show the modulus of the aver- 
age transverse velocities in a few bins of the modulus of the longitu- 
dinal velocity. We see that for TreePM, the magnitude of transverse 



motions drops rapidly as we increase the force softening length. 
We also note that for the adaptive TreePM, the magnitude of trans- 
verse motions is similar to that seen with the TreePM when e for 
the TreePM coincides with the emax for the ATreePM. 

Wc conclude the discussion of this test by noting that our re- 
quirement of e > fij, where fij is the local inter-particle separa- 
tion, ensures coUisionless behavior in evolution of the system. 

5.2 Convergence with n„ and relevance of the Ve term 

For the following discussions we run a power-law model with in- 
dex n = —1.0. The spectrum is normalized at a an epoch when 
the scale of non-linearity r„; = 6.0 in grid units. The scale of non- 
linearity is defined as the scale at which the linearly extrapolated 
mass variance, defined using a top hat filter, is unity a{rni,z) = 1. 
The simulations are done with 128^ particles and a 128^ grid. We 
choose Vs — 1 grid length. Other simulation parameters are as de- 
scribed in the subsection §3.1 on the TreePM code. We assume that 
the Einstein-de Sitter model describes the background universe. In 
this case self-similar evolution of quantities for power law models 
provides an additional test for simulation results. 

5.2.1 Clustering Properties 

We compute the volume averaged 2— point correlation function ^. 

i{r) = ^ r iix)x'dx (26) 
I" Jo 

5 is the two point correlation function (Peebles 1980). To compute 
^ from simulation output, we take 15 independent random subsets 
of 10"' particles each. We then estimate the average value of f over 
these subsets. The maximum and the minimum values of ^ in these 
subsets are our estimate of the errors. Figure 5 shows at an epoch 
when r„i = 6.0. TreePM (left column) with e = 0.025, 0.05, 0.1, 
0.2, 0.5 are shown in black, red, blue, ochre and magenta respec- 
tively. ATreePM with Ve (middle column) and ATreePM without 
Ve (right column) with n„ = 16, 32, 48 are shown in blue, red, 
black lines respectively. The lower row is a zoom-in of ^ so as to 
highlight the differences at small scales, due to the variation in e 
and n„ in different runs. 

One can show that we require a minimum of 12 neighbors to 
solve for Eqn.(20). The Ve term is important in over-dense regions 
where e^ <^ emax, and is a rapidly varying function of density, one 
therefore requires a reasonably large n„ to compute it accurately. 
Use of a small Un leads to a noisy estimate of Ve and Figure 5 
illustrates this point® . With n„ = 16, f deviates significantly from 
our expectations. Disagreement is worse for ATreePM with Ve, as 
this is the case where a poor estimate of the extra term leads to 
larger errors. If we use a larger n„, we expect to see convergence 
at some point. Curves for n„ = 32 and 48 agree with each other 
(within error-bars) for ATreePM with Ve indicating that evalua- 
tion of the extra term is stable. It also indicates that this extra term 
compensates for the use of a larger number of particles, otherwise 
typical softening length is expected to be larger for higher rin and 
this should lead to lowering of the clustering amplitude at small 
scales. This is very clearly illustrated for TreePM (left column in 
Figure 5) where increasing e reduces ^ monotonically. We also see 
this effect for ATreePM without the Ve term where the amplitude 
of clustering is much smaller for n„ = 48 as compared to n„ = 32 

® Also see Figure 2. 




Figure 5. Top Row: This figure plots the volume averaged 2-point correlation function § as a function of scale at r„j = 6.0 for TreePM (left) and ATreePM 
with Ve (middle) and ATreePM without Ve (right). The four curves (black, red, blue, ochre, mauve) for TreePM are for runs with softening lengths of 
e = Jo'io'fo'^''^' "'^ ATreePM runs the three curves (black, red, blue) are for n,i = 16, 32, 48. Bottom Row: zoomed in plots of Top Row. Grey 
surface in the ATreePM plots is § bounded by error-bars for TreePM with e = 
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Figure 6. This figure tests self similar evolution of the scaled 2-point con'elation function ^{r/r„i) by plotting it at the epochs when r„i = 2.0, 4.0, 6.0 
(black, red, blue lines)for TreePM (left) ATreePM with Ve (middle) , ATreePM without Ve 



and the difference between the two is more than twice that seen for 
ATreePM with the Ve term. Such a behavior is expected, and was 
also seen by Price & Monaghan (2007). They noted that the in- 
crease in force softening length with n„ leads to a bias in the force 
at small scales. The Ve term corrects for this and biasing of force is 
less important. Since the distribution of particles gets more strongly 
clustered with time, we expect over- softening of the force field to 
degrade further evolution for the TreePM and the ATreePM without 
Ve, beyond a certain epoch. 

We have plotted the scaled 2-point volume averaged correla- 
tion ^{r/r„i) at small scales in Figure 6 for TreePM with e — 
0.025 (left panel), ATreePM with the Ve term (middle panel) and 
ATreePM without the Ve term (right panel). We used simulations 



with n„ ~ 32 for both the ATreePM runs shown here. ^ has been 
computed at epochs when the scale of non-linearity rni ~ 2, 4 and 
6 (black, red and blue lines respectively). Since the only scale in 
power-law models is the scale of non-linearity r„i, introduced by 
gravity, one expects ^ to evolve in a self-similar manner. The scale 
of self-similarity, r^s, for a given epoch is is determined by finding 
the scale at which ^ matches with the ^ of the earlier epoch within 
the error bars computed in the manner described above. 

We illustrate numerical convergence properties of the three 
codes in table 1. Here, we list the scale Vcon and Vss at different 
epochs. For ATreePM Vcon is the scale beyond which ^ with the 
given n„ matches with a reference ^ with n„ — 48. For TreePM 
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Table 1. The following table shows the variation of Vcon and r^s with time for TreePM and the two variants of ATreePM as a function of softening 



Runs 


Tcon 

Tnl = 6.0 


Tcon 

Tnl = 4.0 


rnl = 2.0 


Vss 

r„i = 6.0 


Vss 

rnl = 4.0 


TreePM (c = ^) 


1.2 


().% 


0.61 


I.I 


0.67 


TreePM (£ = ^) 


0.61 


0.44 


0.36 


0.42 


0.15 


TreePM (e = yg) 


0.25 


0.20 


0.10 


0.18 


0.16 


TreePM (e = §§) 


0.10 


0.10 


0.10 


0.10 


0.17 


TreePM (e= ||) 








0.10 


0.18 


ATreePM (Ve,n„ = 32) 


0.10 


0.15 


0.30 


0.27 


0.31 


ATreePM (Ve,n„ = 48) 








0.42 


0.47 


ATreePM (n„ = 32) 


0.27 


0.28 


0.38 


0.54 


0.42 


ATreePM (n„ = 48) 








0.59 


0.50 



this is the scale where ^ with a given force softening length matches 
with the reference ^ with e = 0.025. 

We find that both variants of ATreePM converge with time, 
the convergence being more rapid for ATreePM with Ve as com- 
pared to ATreePM without Ve. The scale of convergence rcon is 
also smaller for ATreePM with Ve. On the other hand self-similar 
behavior is degraded with evolution for ATreePM without Ve as 
rss increases with time, whereas for ATreePM with the Ve term 
the evolution is self-similar over a larger range of scales. 

Thus we may conclude that the ATreePM with the Ve has well 
defined numerical convergence as we vary n„, and the evolution of 
power law models for this code is self-similar over a wide range of 
scales. This sets it apart from the ATreePM without the Ve term, 
where the convergence with rin is not well defined and the evolu- 
tion of a power law model is self-similar over a smaller range of 
scales. Better match with the fixed resolution TreePM code is an 
added positive feature for the code with the Ve term. This raises 
obvious questions about AMR codes, where no such term is taken 
into account as we increase the resolution of the code. We comment 
on this issue in the Discussion section. 

We briefly comment on the convergence and self-similar evo- 
lution in the fixed resolution TreePM code. As seen in Table 1, we 
find that at early times the scale above which we have self-similar 
evolution is almost the same for e < 0.2. This may indicate dis- 
creteness noise, or it may be due to two body collisions during early 
collapse (Splinter et al. 1998; Melott et al. 1997; Joyce, Marcos, 
& Baertschiger 2008; Romeo et al. 2008). At late times, the scale 
above which evolution is sclf-siinilar is r.^., ^ 2e. This indicates 
that any transients introduced at early times by discreteness, etc. 
have been washed out by the transfer of power from large scales to 
small scales (Little, Weinberg, & Park 1991; Evrard & Crone 1992; 
Bagla & Padmanabhan 1997; Bagla & Prasad 2008). Thus one gets 
an iinproved self-similar evolution at late times in the fixed resolu- 
tion TreePM with the use of a smaller softening length. 

On the other hand we see that for TreePM, rcon increases with 
time. As defined above, this is the scale at which ^ obtained with a 
given value of e matches with the value obtained using e = 0.025. 
We find that at the last epoch rcon > 2e, whereas at early times 
the limit is rcon > e. This trend is contrary to that seen with the 



ATreePM for which rcon decreases with time, the rate at which it 
comes down being faster for the ATreePM with the Ve. We would 
like to point out that we do not probe ^ for scales r < 0.1 as the 
number of pairs with a smaller separation is often too small to get 
reliable estimates of ^. 

Another remarkable fact is that for the TreePM code, rss < 
rcon - Thus we have the desired behavior in terms of evolution even 
though we are yet to achieve numerical convergence at the relevant 
scales. Such a problem does not exist for the ATreePM code. 

Lastly, as we show below, the ATreePM code is very good at 
resolving highly over-dense regions very well. This would appear 
to be in contradiction with the slightly lower two point correlation 
function. The reason for the slightly lower correlation function is 
that the ATreePM code does not resolve small haloes, in particu- 
lar those with fewer than Un particles (see the discussion of mass 
function of collapsed haloes). Indeed, the number density of small 
haloes is severely under-estimated in the ATreePM simulations and 
we believe that this is the main reason for a weaker two point cor- 
relation function. 



5.2.2 Collapsed Haloes 

We now look at another global indicator, namely the mass function 
of collapsed halos, to compare the performance of the TreePM and 
the ATreePM codes. The number density of halos N{M)dM in the 
mass range (Af, M + dM) is plotted in Figure 7. Halos were iden- 
tified using the Friends-Of-Friends algorithm with a linking length 
/; =0.1 in grid units and only haloes with at least 8 particles were 
considered. 

For TreePM we see that as we increase the force softening 
length there is a clear change in the mass function. The mass func- 
tion converges for all masses with e < 0.2. However, the mass 
function for e — 0.5 shows fewer haloes of up to about 10^ par- 
ticles. Thus the effect of a large softening length is seen in mass 
function up to fairly large mass scales. This qualitative behavior 
is insensitive to the choice of linking length, though the scale of 
convergence is smaller if we choose a larger linking length. 

The ATreePM does not sample the haloes with fewer than n„ 
members properly. This is expected as such over-densities are not 



12 Bagla and Khandai 




Figure 7. Mass function for TreePM (left), ATreePM with Ve (middle) and ATreePM without Ve (right). For TreePM t = i Yq i are drawn in black, 
red, blue lines respectively. ri„ = 32, 48 are drawn in black and red lines for ATreePM. The aiTows represent rin in mass units. The black dot-dashed curve 
on every plot is the Press-Schecter mass function, shifted vertically by multiplying by a factor of 1.5. The black dashed line in ATreePM is for TreePM with 




Figure 8. This figure shows slice from simulation where the most massive halo resides (see top left comer). Top panel is for TreePM with e = (left) and 
ATreePM with e = Bottom panel is for ATreePM with Ve (left) and ATreePM without Ve (right). n„ = 32 was taken for ATreePM. 
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Figure 9. This figure is a zoom in of Fig. 8 and shows the region of simulation where the most massive halo resides (see top left comer of Fig. 8). Top panel 
is for TreePM with e = (left) and ATreePM with e = ^. Bottom panel is for ATreePM with Ve (left) and ATreePM without Ve (right). n„ = 32 was 
taJcen for ATreePM. 



sufficient to reduce the softening length from tmax- At larger scales 
we have a convergence for n„ = 32 and n„ = 48, where the 
mass function agrees with that for TreePM with e = 0.025. Mass 
function for n„ = 16 deviates significantly from the other two 
for ATreePM, deviation being larger for the ATreePM with the Ve 
term. We beUeve that this is due to the noisy estimation of the force 
softening length and the Ve term, and resulting errors in force. 

As a reference the Press-Schecter mass function is also plotted 
in all three panels (black dot-dashed line). However it has been 
shifted up vertically by multiplying N{M)dM by a factor of 1.5 
for ease of comparison. 

Distribution of particles in a simulation is a useful represen- 
tation for comparing large scale features. The distribution of par- 
ticles in a thin slice is shown in Figure 8. We have shown the 
distribution for four simulations: TreePM with e = 0.025 (top- 
left), TreePM with e = 0.5 (top-right), ATreePM with the Ve term 
and Un = 32 (lower-left) and ATreePM without the Ve term and 
n„ = 32 (lower-right). The distribution of particles is for the last 



epoch, corresponding to r„i — 6 grid lengths. This slice contains 
the largest halo in our simulation on the top left comer. We note 
that the large scale distribution is the same in all cases, indicating 
that the Adaptive TreePM is not changing anything at large scales. 

We zoom in on the region around the largest halo in Figure 9, 
where the distribution of particles is shown from the set of simula- 
tions used in Figure 8. Distribution of mass at scales larger than a 
gid length is the same in all simulations but at small scales we be- 
gin to sec some differences between the distribution of particles in 
different simulations. TreePM with e = 0.5 differs most from the 
other three, in that it does not resolve small scale structures. This 
happens due to the relatively large force softening length that in- 
hibits collapse if the expected size of the collapsed halo is smaller 
than e. The notable differences between the TreePM (e = 0.025) 
and the ATreePM slices are as follows: 

• Small haloes in the region away from the large haloes are more 
compact for TreePM than for ATreePM. This is perhaps caused 
by the ATreePM not having a constant e and it is likely that in 
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Figure 10. This figure shows the largest halo at r^i = 6.0 for TreePM with e = ^ (left column), ATreePM with Ve and n-n = 32 (middle column) and 
ATreePM without Ve and rin = 32 (right column). The first row contains all the particles in the halo. The second, third and forth rows represent particles in 
the halo which have a density greater than 5 x 10^, 10^ and 2 x 10^ respectively. 
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Table 2. This table shows the number of particles ripart above a density 
threshold (as in Figure 10) retained in the core of the most massive halo 

(O 



Runs 


"^part 

p > o X K)-^ 


l> > 10 " 


p 2 X lO"' 


TreePM (e = ^) 


320 


106 


18 


ATreePM (Ve, n„ = 32) 


568 


297 


125 


ATreePM (n„ = 32) 


374 


169 


18 



these clumps the value of e is larger than 0.025 that is used in the 
TreePM. 

• Location of sub-structure in the central halo differs somewhat 
between the different simulations. 

Figure 10 shows the inner regions of the halo seen in Figure 9. 
We show the inner parts of the halo in the four cases (TreePM (e = 
0.025), ATreePM with the Ve term (n„ = 32), ATreePM without 
the Ve term (rin = 32). The top row shows all the particles in 
the inner parts of the halo. The second row shows all the particles 
with an SPH over-density estimate (computed with n„ — 32 and 
using the kernel specified in Eqn. 5) of greater than 5 x 10*, the 
third row shows the particles with an over-density of greater than 
10^, and, the last row shows particles with an overdensity larger 
than 2 X 10^. We see that the ATreePM with the Ve term shows 
the most pronounced central part of the halo in the top panel, and 
this impression gains strength as we move to lower panels. Table 2 
gives us the number of particles in the various panels of Fig. 10. 
TreePM with e = 0.025 manages to trace some of the highly over- 
dense substructure seen in the ATreePM simulations, though with 
a much sinaller number of particles in these clumps. The ATreePM 
with the Ve term does better than the one without, particularly at 
the highest over-densities used here. Indeed ATreePM with the Ve 
term retains nearly 7 times as many particles in the halo core as 
compared to TreePM and ATreePM without the Ve term. The same 
thing is also reflected in Figure 5 where the ATreePM without the 
Ve term is seen to underestimate clustering at small scales. It is 
noteworthy that the size of highly over-dense structures in the core 
of this halo are much bigger than the force softening length for the 
TreePM. 

It is notable that the ATreePM is able to resolve highly over- 
dense regions while taking much less time than the fixed resolution 
TreePM. Thus we have the added performance at the cost of fewer 
resources. 



5.2.3 Dynamics within Collapsed Haloes 

The equation of motion for the adaptive code has the additional Ve 
term, and it is important to test whether this term leads to a change 
in dynamics in highly overdense regions. We test this by studying 
the dynamics in central cores of the largest haloes in the simulation. 
We study the ratio {2T)/{U) at a scale corresponding to r2oo/5, 
with the averaging done arotmd the centre of mass of the halo. We 
chose the reference scale in this manner as the density profiles for 
haloes with different softening length differ considerably and using 
a high density contour to define the radius leads to very different 
physical scales. Figure 1 1 shows {2T) / (f/) as a function of soften- 
ing length e for the fixed resolution TreePM simulations. This has 
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Figure 11. Ratio of two times the kinetic energy and the potential energy, 
2T/U, for the 5 largest halos as a function of softening length for TreePM 
(black hnes). Open squares and filled triangles are for ATreePM with and 
without Ve respectively. 

been done for the five largest haloes in the simulation. We see that 
as we go to a larger softening length, the ratio becomes larger than 
unity, indeed for e = 0.5 the ratio is closer to two for one of the 
haloes. This is likely to happen if the size of the over-dense core 
is comparable to the softening length, and indeed this is the case 
for some of the haloes used for this plot. The same plot shows the 
ratio as seen in the adaptive code (with the Ve term) with empty 
squares, and as filled triangles for the adaptive code without the Ve 
term. These circles are plotted at small values of e, not used for 
TreePM runs. The values of e used for positioning these symbols in 
the plot has no relevance to the simulation, and this has been done 
purely for the purpose of plotting the values. We see that the val- 
ues for the ratio in both the adaptive codes correspond closely to 
the values of the ratio seen in fixed resolution TreePM simulations 
with a small e. Thus we may conclude that the additional term in 
the Adaptive code is not leading to any significant changes in the 
dynamics, as compared to the fixed resolution codes. 



6 DISCUSSION 

In this paper we have introduced the first cosmological N-Body 
code that has a continuously varying force resolution. This is based 
on a well defined formalism that ensures energy and momentum 
conservation. An important aspect of this formalism is that it mod- 
ifies the equation of motion at scales below the force softening 
length. Our implementation of the formalism uses methods devel- 
oped for SPH codes, this is required as we need to assign quantities 
like density to particles. We have described our implementation of 
the code in detail here. We have given a brief summary of errors 
in force for the Adaptive TreePM code with the new parameter n„. 
Given that this code is based on the TreePM code, it inherits the 
errors in force with respect to all the other parameters. As we have 
seen in earlier work (Bagla 2002; Bagla & Ray 2003; Khandai & 
Bagla 2008), errors in force can be minimised to a fairly low level 
for the TreePM with a judicious choice of parameters. 

We find that the time taken by the Adaptive TreePM (for n„ = 
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32) is comparable to that taken by the fixed resolution TreePM code 
if the latter uses a force softening length that is about 1/20 of the 
mean inter-particle separation. 

One of the key reasons for considering adaptive resolution is 
to avoid two body collisions. Cosmological simulations require col- 
lisionless evolution. We test the Adaptive TreePM by simulating 
collapse of an obUque plane-wave. We find that unhke the fixed 
resolution TreePM code that leads to large transverse motions, a 
clear sign of two body collisions, the Adaptive TreePM code does 
not have any significant transverse motions. (See Figures 3, 4) 

We have compared the performance of the Adaptive TreePM 
with that of the fixed resolution TreePM code for power law initial 
conditions. We used a power law model and the Einstein-de Sitter 
background as it allows us to test codes by requiring self-similar 
evolution of the correlation function. We find that for n„ = 32, the 
resulting density distribution in the Adaptive TreePM is compara- 
ble to that seen in the highest resolution TreePM in many respects. 
The correlation function for the two match at all scales and the mass 
function of collapsed haloes matches for haloes with more than n„ 
particles. The Adaptive TreePM with the Ve term in the equation of 
motion does much better than the TreePM as well as the Adaptive 
TreePM without the extra term in resolving highly over-dense cores 
of massive haloes. It is noteworthy that the Adaptive code takes less 
time than the TreePM it has been compared with (see fig. 1). 

We have tested the codes by looking for convergence in results 
as we modify the key parameters that describe the force resolution. 
We find that the fixed resolution TreePM code converges slower 
than expected from, for example, self-similar evolution of cluster- 
ing. The Adaptive TreePM code with the Ve term converges fairly 
quickly and offers an effective dynamical range that is sUghtly 
smaller than the fixed resolution TreePM that takes almost the same 
time to run. 

Given our analysis of errors, and the comparison of the per- 
formance of the ATreePM code with and without the extra term in 
the equation of motion, the natural conclusion is that we require the 
extra term in order to obtain low errors and numerical convergence. 
This raises the obvious question about the AMR codes where no 
such extra term is used. Further, in most AMR codes, resolution 
is increased when there are of order 10 particles in the lower res- 
olution elements. We find that the errors are minimized when this 
number is around 20 even when the extra term is not taken into ac- 
count. It is not clear how serious these issues are, given that AMR 
codes have been developed and tested in a variety of ways over the 
last three decades, but it is a concern®. 

We have established in this paper that an adaptive resolution 
code for evolution of perturbations in collisionless dark matter can 
give reliable results for a range of indicators from clustering prop- 
erties to mass function of collapsed haloes and even get the internal 
dynamics of collapsed haloes tight. Further, we show that such a 
code is efficient in that it is faster than fixed resolution codes that 
give us comparable force resolution in highly over dense regions. 
Such a code is very useful as it allows us to probe clustering at small 
scales in a reliable manner. Studies of box-size effects have shown 
that large simulation boxes are required in order to limit the ef- 
fect of perturbations at larger scales that are not taken into account 

® One aspect where AMR codes do relatively poorly is in getting the halo 
mass function at the low mass end (O'Shea et al. 2005; Heitmann et al. 
2008). As one can see in figures in the papers cited here, the shortfall is 
often spread over a decade in mass of haloes. We do not see any gradual 
decline in the number density of haloes in ATreePM up to the cutoff set by 
n„. 



(Bagla & Prasad 2006; Bagla, Prasad, & Khandai 2008). In such 
a situation an adaptive code provides us with a reasonable range 
of scales over which the results can be trusted. We expect to use 
this code to address several issues related to gravitational cluster- 
ing in an expanding universe. Wc also plan to revisit issues related 
to density profiles of collapsed haloes. 

Given the conclusions listed above, we feel that the Adaptive 
TreePM methods represents an exciting development where we can 
set aside worries about the impact of coUisionaUty. The relative 
speed of the adaptive code also makes this a more pragmatic option. 
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APPENDIX A: EXPRESSIONS FOR CUBIC SPLINE SOFTENED POTENTIAL AND FORCE KERNELS 



One can integrate the Poisson equation for a given kernel to obtain the softened two-body potential kernel. For the case of the cubic spline 
kernel we have the expression for the two-body potential kernel 0(u, e) and the regular two-body force f (u, e) = — Vu0 : 



(u,e) 



f(u,e) 



-^[|-!(7)^+(7)T- 



+¥(7) 



(7)1 



< f < 0.5 
iS, 0.5 < 2i < 1.0 

1.0 < f 

< ^ < 0.5 

0.5 < f < 1.0 

1.0 < ^ 



(Al) 



(A2) 



For variable or local smoothing length ti = e(rj) there will be an additional term given in eqs.(10-12) for which one needs to compute 
(dW/de), (dW/du) and (90/9e) For the cubic spline kernel these expressions are: 
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